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Automated Method for identifying landmarks within an image of the brain 
Field of the invention 

The present invention relates to methods for automatically identifying 
landmarks within images of a brain. 

5 Background of Invention 

The Talairach transformation is widely used for analysis of neurological 
images. It involves identifying eight landmarks, which are used to define a co- 
ordinate system. The Tailairach landmarks subdivide the brain into 12 
cuboids, and the Tailairach transformation is to warp the images within each 
cuboid linearly. In this way the brain images are normalised by a three- 
dimensional piece-wise linear warping. This scheme has several applications, 
in particular because it makes it possible to compare neurological images 
from different individuals. One improvement on this scheme, while following 
its conceptual rationale, is the improvement of the definitions of the 
landmarks, to give "modified Talairach landmarks" (as defined in the article 
"Modified Talairach Landmarks", W. L. Nowinski, Acta Neurochirurgica, 2001 , 
143, p1 045-1 057, the disclosure of which is incorporated herein by reference). 

20 A principal advantage of the Talairach transform is its simplicity. Although 
numerous non-linear image registration techniques are known, in principle 
providing a higher accuracy, the non-linear techniques have limitations which 
make them difficult to use beyond a research environment. In particular, a 
prohibitively high computational time is required. Whereas the Talairach 

25 transformation can be performed in less than a second in a standard personal 
computer, some non-linear methods require days of computation. More 
fundamentally, the non-linear methods are conceptually complex, and must 
be treated as "black boxes", which limits their clinical acceptance. 
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Furthermore, there is no established methodology for validation of the non- 
linear registration techniques. 

One drawback of the Talairach transformation, however, is that the 
5 identification of the landmarks has not so far been automated reliably, so that 
when using conventional software which employs the Talairach transformation 
the landmarks still have to be identified by user interaction (e.g. R. W. Cox, 
"AFNI: Software for, Analsysis and Visualization of Functional Magnetic 
Resonance Neuroimages", Computer and Biomedical Research, 1996, 29, 
10 p1 62-1 73). Quite apart from the time the interactive identification of the 
landmarks takes, different individuals are liable to locate the landmarks in 
slightly different positions, which reduces the robustness of the method. 

Summary of the Invention 

15 The present invention aims to provide new and useful techniques and 
apparatus for identifying automatically landmarks within neurological images, 
and in particular the Talairach landmarks or modified Talairach landmarks. 

The invention makes use of the concepts "foreground voxels" and 
20 "background voxels". "Binarising" a set of voxels means to divide the voxels 
into these two categories. 

Furthermore, certain voxels are "object voxels". Object voxels are voxels in a 
location which corresponds to a physical identity (cerebral structure), such as, 
25 in the context of this invention, the union of grey matter and white matter. 

In general terms, the invention proposes that the location of one or more 
landmarks is obtained by the following steps: 



3 

(a) a region of interest in a plane within an MR image of a brain and 
containing the landmark(s), is identified; 

(b) the region of interest is binarised into foreground and background 
voxels based on at least one threshold selected using anatomical knowledge, 

5 (c) a set of object voxels is identified from the foreground voxels, 

excluding voxels which were only classified as foreground voxels due to 
proximity of cortical and non-cortical structures, 

(d) object voxels are identified from the background voxels due to the 
partial volume effect and morphological erosion/opening, and 

10 (e) an automatic process is then carried out to identify one or more 

landmarks in the modified binarised image. 

The step of reclassifying voxels may be performed in two stages: (i) a first 
stage of morphological processing, and (ii) a second step of restoring voxels 
which have been incorrectly reclassified in the morphological processing. 

15 Anatomical knowledge may also be used in the reclassifying step, e.g. by 
using the expected shapes of cortical and/or non-cortical structures to modify 
the classification. 

The threshold is preferably selected by the steps of: 

(i) using prior knowledge about the image to derive an intensity range of 
20 voxels in the said region of interest; 

(ii) obtaining a frequency distribution of the intensities within the said intensity 
range of voxels within the said region of interest; and 

(iii) using the said frequency distribution to derive the threshold. 
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The intensity threshold may be selected by minimising a function which is a 
sum of the variances of the intensities below and above the threshold. 

Optionally, this function may be a weighted sum defined based on two 
constants W 1 and W 2 . This is referred to here as a "range constrained 
5 weighted variance method". 

For example, labelling the possible values of voxel intensity by an integer 
index i and their respective frequencies by h(i) t and writing the lower and 
upper intensities respectively as n ow and r high , the weighted sum is given by 

o*clwv<Wx>wj = maxW^F, + p r (c 2 )D(c 2 )^ 2 ), 

10 where Pr(.) denotes the class probability (PrCC,) = ]>]/*(/) and 

i=r tow 

r hi g h 

Pr{C 2 ) = 2^(0 ).. and D(d) and D(C 2 ) are given by: 

'high 

D{C X ) = G" 0 - n T ) 2 and D(C 2 ) = (//,- ju T Y » where v T = £ f x *(0 ■ 

>"o = S^^Oand //, = ]T/x/i(0. 

The steps (a) to (d) may be performed repeatedly for different landmarks. One 
15 or more landmarks may be located for each of the starting planes of the MR 
image. 

In the case of identifying the I landmark (which is a position on a coronal 
plane) the left and right halves of the brain may be treated partly separately, 
such that the I landmark is identified using the half of the brain which appears 
20 to have been segmented more accurately. 



The invention may be expressed as a computer-implemented method, or 
alternatively as a computer system arranged to perform such a method. 
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Brief Description of The Figures 

Preferred features of the invention will now be described, for the sake of 
illustration only, with reference to the following figures in which: 

Fig. 1 is a flow chart of a method which is an embodiment of the 
5 invention; 

Fig. 2 is a flow chart showing in more detail sub-steps of step 3 of the 
flow chart of Fig. 1; . 

Fig. 3, which is composed of Figs. 3(a) to 3(d), shows an example of 
performing the steps to identify the A, P, L and R landmarks in the method of 
10 Fig. 1; 

Fig. 4, which is composed of Figs. 4(a) to 4(d) , shows an example of 
performing the steps to identify the S landmark in the method of Fig. 1; and 

Fig. 5, which is composed of Figs. 5(a) to 5(d), shows an example of 
performing the steps to identify the I landmark in the method of Fig. 1 

15 

Detailed Description of the embodiments 

Referring firstly to Fig. 1, the steps of a method which is an embodiment of the 
invention are shown. 

20 

In step 1, a dataset which is a neuroimage (i.e. an image of a brain) is loaded. 
This image is a three dimensional data, typically obtained from a number of 
scans in different respective planes. 

25 From this data, the midsagittal plane (MSP) is determined. This is preferably 
done using the method disclosed in WO02/069827, "Method and apparatus 
for determining symmetry in 2D and 3D images", by Hu and Nowinski (the 
disclosure of which is incorporated herein by reference). However, the 
invention is not limited in this respect, and any other technique for determining 
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the MSP may also be applied. Indeed, it would also be possible within the 
scope of the invention for the input data to specify the MSP. 

The coordinates of the anterior commissure (AC) and posterior commissure 
5 (PC) are then determined automatically. This can be done by the method 
disclosed in WO02/43003, "Methods and apparatus for processing medical 
images", by Nowinski and Thirunavuukarasuu (the disclosure of which is 
incorporated herein by reference), although once more the invention is not 
limited in this respect. 



In step 2 of Fig. 1, the data is normalised to occupy a predefined gray-scale 
range. According to the standard radiological convention, we write the 
coordinate system as (x,y,z), where the x-axis is from the subject's right to left, 
the y-axis is from anterior to posterior, and z from superior to inferior. Thus, an 

15 xz-plane is a coronal plane, a yz-plane is a sagittal plane, and an xy-plane is 
an axial plane. Let g(x,y,z) denote the gray level of the input data at a voxel at 
position (x,y,zj, and let g 0 and gr* denote the grey levels such that one percent 
of the voxels have an intensity less than go and one percent of the voxels 
have an intensity greater than g?. Then, we obtain a normalised gray level 

20 g(x,y,z) which, for a given position (x,y,z) is given by: 
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g(x 9 y,z)=0 



if g(*,y,z)<g 0 



g(x 9 y 9 z) = 



g(x 9 y 9 z)-g 0 
S\ ~~ So 



if go <g(*>y>z)<gi 



g(x,y,z) =255 



'f g(x 9 y,z)>g x 



25 



Each slice of the normalised data has its own co-ordinate system (u,v) where 
u is the horizontal direction and v is the vertical direction. 
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In step 3, the position of the A, P, L and R landmarks is located. This is done 
by the series of steps shown in Fig. 2. 

Firstly, in step 3.1 a region of interest (ROI) is defined. These may be the 
voxels within the skull. These voxels can be determined by the following 
5 steps: 

(a) A histogram-based thresholding method is used to binarise the AP plane 
(as used for example in M. E. Brummer, R. M. Mersereau, R. L. Eisner, R.J. 
Lewine, "Automatic detection of brain contours in MRI data sets", IEEE 
Transactions on Medical Imaging 1993; 12(2), p153-166). 

10 (b) A morphological closing operation is performed using a 3x3 structuring 
element (SE) to connect small gaps in the ROI. 

(c) The largest connected component is identified, and the holes within the 
component are filled. 

Fig. 3(a) shows the AP plane in a typical example of the use of this method. 
15 Fig. 3(b) shows the corresponding ROI. 

In step 3.2 an optimum threshold is determined, based on the range- 
constrained weighted variance thresholding method. This includes following 
steps, which are explained in a separate patent application by two of the 
present inventors: "Methods and apparatus for binarizing images", Singapore 
20 patent application 200307531-4, by Q. M. Hu, Z. Hou, and W.L. Nowinski, 
which was still unpublished at the priority data of the present application, and 
the disclosure of which is incorporated by reference. 

Firstly, prior knowledge of the image is used to define an ROI which is a 
subset of the image. This process can be done by whatever means, either 
25 automatic, semi-automatic, or even manual. Then an analysis is performed on 
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the frequency of occurrence of intensities within the ROI, and a range of 
frequencies is defined, again using prior knowledge. 

For example, without losing generality, we denote the image to be processed 
5 as f(x), where f(x) is the gray level at a voxel labelled x. It is further supposed 
that the processed image has L gray levels denoted by r,. where / is an integer 
in the range 0 to L-1 and r Q <r, <...r L _ x . It is also assumed that the object of 
interest has higher intensity values than the background. Suppose that due to 
prior knowledge or test we know that the proportion of the ROI which is 
10 occupied by the object is in the percentage range per 0 to pen. 

Let h(i) denote the frequency of gray level r h and let H(i) denote the 
cumulative frequency which is ^Ki'), where /' is an integer dummy index. 
Considering two values of / written as m and n, the frequency of intensities in 

n 

15 the range r m to r n is . Thus, we can use per 0 to calculate a gray level 

n ow , such that we are sure that all the voxels having lower intensity represent 
background. n ow can be written as: 

r, ow = m / n {' I H (0 * P er o } ■ 0 ) 

20 Similarly, we can use pen to calculate a gray level r high such that we are sure 
that all voxels having higher intensity represent the object: 

r high = min{/ 1 H(i) > per, } . (2) 

Let r k fall within the range n ow to r hi g h , and suppose that the voxels of the ROI 
25 are in two classes Ct and C 2l where C* is voxels of the background class and 
consists of voxels with gray levels n 0 w to and C2 is voxels of the object class 
and is composed of voxels with gray levels r k +1 to r high . The range-constrained 
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weighted variance thresholding method maximises the "weighted between- 
class variance" defined as: 

0*cwv W >W 2 ) = max (Pr(C, )D(C, Wi + Pr(C 2 )D(C 2 )W 2 ), 

where Wi and tV 2 are two positive constants selected by the user and 
5 representing the weights of the two respective class variances, Pr(.) denotes 
the class probability, i.e. 

Pr(C,)= 2>(f), Pr(Q)= £ /»(/). 

and D(Ci) and D(C 2 ) are given by: 

r h-,gh 

^(C 1 ) = C" 0 -// r ) 2 and Z)(C 2 ) = Ct/ I -// r ) 2 , where ju T = ^ixhQ), 

10 /*o = xA(0 and //j = J/x^(/). 

When is bigger than tV 2> background homogeneity is emphasised. 

Step 3.2 may be done by specifying per 0 and pen to be 14% and 28% 
respectively. The optimum threshold is denoted as 0, . 

In step 3.3, we segment the AP plane by assigning voxels to foreground if 
15 they are bigger than 0 X% and otherwise assigning them to background. The 
binarised image is denoted as BWAP1(u,v). 

In steps 3.4 we reclassify the voxels: firstly by a morphological processing, 
then processing using anatomical knowledge, and finally performing a 
restoring step. 

20 The sub-steps of the morphological processing and processing using 
anatomical knowledge are as follows: 



(a) A distance transform of the ROI is performed using the 2-3 metric (in 
this metric the distances between any two voxels is determined by defining 
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the shortest path of voxels between them, and adding the distance increments 
along this path. The distance increments between two voxels which are 
nearest neighbours in a direction parallel to one of the axes (4-connected 
nearest neighbours) is taken as 2, and the distance between two voxels which 
5 are nearest neighbours diagonally (8-connected nearest neighbours) is taken 
as 3) and the distance codes are converted into distance indices by the 
method of Hu Q.M., "Two and three-dimensional skeletonization", WO 
02/058008). The maximum distance index is denoted as maxDSkull. 

(b) A morphological opening operation is performed with a 3x3 SE with 
10 respect to BWAP1(u, v), to obtain BWAP2(u,v). 

(c) A morphological opening operation is performed with a 5x5 SE with 
respect to BWAP1 (u,v) to obtain BWAP3(u,v). 

(d) The foreground components of BWAP3(u,v) are found. For each 
foreground component, its minimum and maximum distance indices are 

15 denoted as minD and maxD respectively. A foreground component is treated 
as an object component when maxD-minD is bigger than a value (e.g. 20) or 
maxD is bigger than a second value (e.g. maxDSkull/2). 

(e) The object voxels are excluded from the foreground voxels of 
BWAP2(u,v). The connected foreground components of BWAP2(u t v) are 

20 found. The number of voxels of each foreground component are denoted as 
nosVoxel. A foreground component of BWAP2(u,v) is categorised as an 
object component only when the shape of the component is not similar to 
meninges. According to anatomical knowledge, meninges have a shape 
similar to the skull and are quite thin. So, when maxD-minD is smaller than 

25 0.1*nosVoxel, the component is highly likely to be a meninges. Otherwise, it is 
classified as an object component. 

The lost object voxels are restored by the following steps: 
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(a) Object voxels far from the skull lost due to the morphological opening 
operation are restored. This is achieved by checking the non-object voxels 
with a distance index greater than maxDSkull/4. If their gray level is bigger 
than 9 X , the voxels are reclassified as object voxels. 

5 (b) Object voxels around the boundaries lost due to the morphological 
opening operation are restored. For each object boundary voxel (an object 
boundary voxel is an object voxel having at least one non-object voxel as one 
of its 8 immediate neighbours), each of its 8 immediate neighbours is 
reclassified as an object voxel if its grey level is greater than 9 X . 

10 Note that this restoration is not performed around the most anterior (i.e. 
minimum v) and most posterior (i.e. maximum v) parts of the straight line 
connecting AC and PC, to avoid the sinus and meninges. Specifically, 
suppose the maximum and minimum v coordinates of object voxels are 
minVap and maxVap respectively, and denote the coordinates of AC on the 

15 AP plane as (acUap,acVap) and the coordinates of the PC on the AP plane as 
(pcUap, pcVap). The restoration is not carried out in the following two regions: 

|u-acUap|<10mm and |v-minVap|<3mm. (3) 

|u-pcVap|<10mm and |v-maxVap|<3mm, (4) 

where |x| stands for the absolute value of x. 

20 (c) Object voxels (that is, all voxels which are physically either GM or WM) lost 
due to the partial volume effect are restored. Since the statistics of GM, WM, 
CSF, air, meninges and bones are not available, the partial volume effect is 
alleviated by reducing the threshold by 10. That is, for any object boundary 
voxel, each of its immediate 8 neighbours is re-classified as an object voxel if 

25 its gray level is bigger than (0,-10). The restoration is not carried out in the 
two regions defined by formulae (3) and (4). 
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In step 3.5, the minimum and maximum v coordinates of the object voxels are 
taken respectively as the v coordinates of the A and P landmarks respectively. 
Similarly, the minimum and maximum u coordinates of the object voxels are 
taken as the u coordinates of the L and R landmarks respectively. Note that 
5 the u-coordinate in the AP plane corresponds to the x-coordinate of the three- 
dimensional volume, and the v-coordinate in the AP plane corresponds to the 
y-coordinate in the three-dimensional volume. 

Figs. 3(c) and 3(d) show the segmented AP plane, and the 4 landmarks 
overlaid over it. The white horizontal lines show the v coordinates of the A and 
10 P landmarks, while the white vertical lines shown the u coordinates of the L 
and R landmarks on the AP plane. 

In step 4 of Fig. 1, we determine the position of the S landmark. The VPC 
plane is a coronal slice perpendicular to both the MSP and the AP plane, and 
it passes through the PC. To determine the position of the S landmark, we 
15 only need to determine its v co-ordinate. The v coordinate of the S landmark 
is the smallest v coordinate of all the cortical voxels on the VPC plane. The S 
landmark is localized by segmentation of a virtual slice aVPC(u,v) with the 
close skull. 

The VPC plane is denoted as VPC(u.v), the coordinates of the PC within 
20 VPC(u,v) are denoted as (pcLI, pcV). 

In step 4.1, a virtual plane aVPC(u.v) is constructed in the following way: 

(a) aVPC(u,v) = VPC(u.v) when v is not bigger than pcV. 

(b) aVPC(u,v) =VPC(u,2pcV-v) when v is bigger than pcV and smaller than 
2pcV. 



25 Figs. 4(a) and 4(b) show a VPC and the corresponding virtual slice aVPC. 
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In step 4.2, the ROI corresponding to aVPC is found. This procedure is done 
in the same way in which the ROI for the AP plane was found above (step 4.1 
above). 

In step 4.3, the optimum threshold 0 2 is determined by the range-constrained 
5 weighted thresholding method, by specifying per 0 and pen to be 20% and 
40% respectively. 

In step 4.4, the aVPC plane is segmented using the optimum threshold 0 2 , by 
the same set of sub-steps explained in step 4.3 above. 

In step 4.5, the threshold 0 2 is adjusted using anatomical knowledge. Since 
10 the vertical line passing through the PC in the vicinity of the S landmark 
should be interhemispheric fissure voxels, 0 2 should be higher than the gray 
levels of voxels on the vertical line segment starting from 2mm above and 
ending 2mm below the object voxels with the minimum v coordinate found in 
step 4.4. If <9 2 is indeed bigger than this value, then it is not changed. 
15 However, if it is lower, it is modified to be 5 plus the maximum gray level of 
the line segment, and the aVPC(u,v) is re-segmented with the modified 
threshold 9 2 , by the same sub-steps as those used to segment the AP plane 
in step 4.3 above, including the same morpological opening operations. 

In step 4.6, lost object voxels are restored. This is done by the following steps: 

20 (a) Object voxels around the object boundaries lost due to the 
morphological opening operation are restored. This is done by, reclassifying 
an non-object voxel which has an object boundary voxel as a nearest 
neighbour and which has a gray level greater than0 2 . 

Note, however, that to avoid the sinus/meninges, the restoration is not carried 
25 out in the region defined by: 
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|u-pcU|<10mm and |v-minVavp|<3mm (5) 

where minVavp denotes the minimum v coordinate of object voxels. 

(b) Object voxels lost due to the partial volume effect are restored. This is 
done by reclassifying an non-object voxel having an object boundary voxel as 
5 a nearest neighbour and which has a gray level higher than (0 2 -1O). The 
restoration is not, however, carried out in the region defined by equation (5). 

The segmented aVPC after restoration sub-steps (a) and (b) is shown in Fig. 
4(c). 

In step 4.7, the v coordinate of the S landmark is the mimimum v coordinate of 
10 all object voxels in aVPC. The v coordinate of S in the aVPC plane is equal to 
its z coordinate in the full three-dimensional volume. 

Fig. 4(c) shows the segmented aVPC slice, and Fig. 4(d) shows the original 
VPC overlaid by a horizontal line indicating the z coordinate (or equivalent^ v 
coordinate) of the S landmark. 

15 In step 5, the position of the I landmark is determined. The VAC plane is a 
coronal slice parallel to the VPC plane, and passes through the AC. Only the 
v coordinate of the I landmark is need, and it is the maximum v coordinate of 
all the cortical voxels on the VAC plane. It is determined by segmenting the 
VAC plane directly constrained by anatomical knowledge. It is assumed that 

20 the maximum difference in the z-coordinate between the AC and the I 
landmark is no more than 50mm. We denote that coordinates of the AC in the 
VAC(u,v) by (acU, acV). 

Specifically, the I landmark is obtained by the following sub-steps: 
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In step 5.1, VAC(u,v) is binarised by assigning all voxels with gray levels 
bigger than 9 2 to be foreground voxels, and the rest as background voxels. 
The binarised image is denoted as BWVAC1(u,v). 

In step 5.2, the region around the AC is connected to make subsequent 
5 "seeding" feasible. Here "seeding" means to find connected components from 
a specified voxel ("seed") with all the voxels in the component being of the 
same type (i.e. background, foreground or object). The region around the AC 
is connected by setting the BWVAC1(u,v) to foreground when |v-acV| is 
smaller than 3mm. 

10 Step 5.3 employs a vertical line passing through the AC to separate the VAC 
into left and right halves. The voxels on the vertical line with a v-coordinate 
greater than acV+3mm are set to background. That is, BWVAC1(u,v) is set to 
background when v is bigger than (acVac+3mm) and u is equal to acU. This 
has the effect that there is foreground separation in the neck region. 

15 In step 5.4, a morphological opening operation is performed on BWVAC1(u,v) 
using a 3x3 SE, to give BWVAC2(u,v). This operation breaks weak 
connections between the cortex and the skull and between the cortex and the 
neck. 

In step 5.5, a morphological erosion operation is performed using a 3x3 SE, to 
20 give BWVAC3(u,v). This operation further breaks the connections between 
the cortex and the skull and between the cortex and the neck. 

In step 5.6, we seed from (acU, acV), to obtain the foreground component. 
Then, we perform a morphological dilation on the seeded foreground 
component with a 3x3 SE, to obtain BWVAC4(u,v). The serial operations 
25 (erosion followed by seeding and dilation) are intended to break strong 
connections between the cortex and the skull, and between the cortex and the 
neck, while retaining the original shape of the cortex. 
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In step 5.7, the maximum value of v for which a foreground BWVAC4(u,v) 
voxel exists having u smaller than acU, is found, and denoted as maxVL. 

In step 5.8, the maximum value of v for which a foreground BWVAC4(u,v) 
voxel exists having u not less than acU, is found, and denoted as maxVR. 

5 In step 5.9, the left half of BWVAC4 (u,v) (i.e. the values of u smaller than 
acU) is restored in two substeps if (maxVL-acV) is smaller than 50mm: 

(a) Firstly, the effects of the morphological opening operation are 
counteracted, by changing any background voxel which has at least a 
foreground voxel of BWVAC4(u,v) as one of its 8 immediate neighbours and 

10 which has a gray level in VAC(u,v) greater than 0 2 to foreground. 

(b) Secondly, the effects of the partial volume effect are counteracted by 
restoring any background voxel which has a foreground boundary voxel (a 
foreground boundary voxel is a foreground voxel, within its 8 immediate 
neighbours there is at least a background voxel) of BWVAC4(u,v) as one of its 

15 8 immediate neighbours and which has a gray level in VAC(u.v) bigger than 
( 0 2 -1 0) to foreground. 

Similarly, the right half of BWVAC4(u,v) (i.e. the half with u not less than acU) 
is restored by two corresponding steps when maxVR-acV is smaller than 

50mm. 

20 Object voxels are all foreground voxels in this case. 

In step 5.10, if both (maxVL-acV) and (maxVR-acV) are smaller than 50mm, 
the v coordinate of the I landmark is obtained as the v component of the 
object voxel in BWVAC4(u,v) with the biggest value of v. If one of (maxVL- 
acV) or (maxVR-acV) is smaller than 50mm and the other is not, the v- 
25 coordinate of the I landmark is obtained as the maximum v coordinate of all 
the object voxels from the side (i.e. the left side or right side) for which the 
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maximum v coordinate of the object voxels is smaller than (50+acV). If both of 
(maxVL-acV) and (maxVR-acV) are bigger than 50mm, the v coordinate of 
the I landmark Is the maximum v coordinate of all the object voxels on the 
side (i.e. left or right) whose maximum object v coordinate is smaller. Note 
5 that the v coordinate of BWVCA(u,v) corresponds to the z coordinate of the 
dataset. 

Fig. 5(a) shows the original VAC, Fig. 5(b) shows the binarised VAC (i.e. 
BWVAC1), Fig. 5(c) shows the processed foreground (BWVAC4), and Fig. 
5(d) shows the v coordinate (or equivalently z coordinate) of the I landmark 
10 overlaid on the original VAC plane. 

Note that the set of sub-steps performed in steps 4 and 5 can be considered 
as corresponding to those shown in Fig. 2 for finding the A, P, L and R 
landmarks. That is, a ROI is identified; a threshold is selected (or, in the case 
of step 5, the threshold used is that same one derived in step 4); a 
15 segmentation is performed; then a reclassification is performed; and finally the 
landmarks are identified. 

Although only a single embodiment of the invention has been described in 
detail, many variations are possible within the scope of the invention as will be 
clear to a skilled reader. 
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